function [ft, fcoi] = interpAt( t, sol, field, indices )
% Interpolate solution
%
global TJ TJsum                             % To get COI
num = size(indices, 1);
nnode = sol.nnode; ft = zeros(size(t,1), num);
nstate = sol.nstate; nG = sol.nstate /2; nctrl = sol.nctrl;
state_e = (nnode+1)*nstate; ctrl_e = state_e + nnode*nctrl;
t0 = sol.time(1:end-1);

if strcmp(field, 'states') == 1
	stateMat = reshape(sol.x(1:state_e),nnode+1,nstate);
    coi = sum(stateMat(:, 1:nG) * diag(TJ), 2) /TJsum;
    for k = 1 : num
		ft(:,k) = lagfun(t0, stateMat(1:end-1,indices(k)), t) * 180 / pi;
    end
    fcoi = lagfun(t0, coi, t) * 180 / pi;   % To get COI
elseif strcmp(field, 'controls') == 1
	ctrlMat = reshape(sol.x(state_e+1 : ctrl_e),nnode,nctrl);
	for k = 1 : num
		ft(:,k) = lagfun(t0, ctrlMat(:,indices(k)), t);
	end
end
%%----------------  lagrange function  ----------------%%
	function fx = lagfun(x0, fx0, x)
		nx0=length(x0); fx = zeros(size(x,1), 1);
		for i=1:nx0
		    p=1;
		    for j=1:nx0
		        if i~=j, p = p.*(x-x0(j))/(x0(i)-x0(j)); end
		    end
		    fx = fx + fx0(i)*p;
		end
	end
end

